%	Compute the canonical example of Section IX of LANCIA, RUSSO, and WORRALL (2024) using the shooting algorithm %
clear all; close all;
tic
warning('off')

cd ..\
cd Matlab_functions

global pi1 pi2 betta delta varxi xaut yaut x_target y_target time 

the=1;

u_function=@(xxx)log(xxx);

gamma = 1;                              % Risk Aversion
pi1 = 0.5;                              % prob state 1
pi2 = 1-pi1;                            % prob state 2
theta = 0;                              % aggregate shock  
e1 = 1+theta;                           % endowment state 1
e2 = 1-(pi1/(pi2))*theta;               % endowment state 2
sigma = 0.1;                            % measure of volatility of endowment
kappa = 3/5;                            % measure of mean of endowment
ey1 = (kappa+sigma*pi2/pi1)*e1;         % endowment state 1 young
eo1 = (1-kappa-sigma*pi2/pi1)*e1;       % endowment state 1 old
ey2 = (kappa+sigma)*e2;                 % endowment state 2 young
eo2 = (1-kappa-sigma)*e2;               % endowment state 2 old
betta = exp(-1/75);                     % individual discount Factor
delta = betta;                          % social planner discount Factor

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Autarky Allocation %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

aaut  = eo1/e1;
baut  = eo2/e2;
xaut=aaut/(1-aaut);  %MRS in state 1
yaut=baut/(1-baut);  %MRS in state 2

varxi=xaut*((1/(1+xaut)))^(((1+betta*pi1)/(betta*pi1)))...
    *((yaut/(1+yaut)))^(((pi2)/(pi1)));   % useful parameter for time-varying allocation

v1aut=(u_function((ey1))+betta*(pi1*u_function(eo1)+(1-pi1)*u_function(eo2))); % intertemporal utility in autarky in state 1
v2aut=(u_function((ey2))+betta*(pi1*u_function(eo1)+(1-pi1)*u_function(eo2))); % intertemporal utility in autarky in state 2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Time-varying allocation %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Longrun of time-varying allocation
tds=@(z)[% FOC state 1
        z(1)-(((varxi.*(1+z(1)).^((1/(betta.*pi1))).*(((1+z(2))./(z(2))))...
        .^(((pi2)/(pi1))))./(1-varxi.*(1+z(1)).^((1/(betta.*pi1)))...
        .*(((1+z(2))./(z(2)))).^(((pi2)/(pi1))))));
         %FOC state 2      
        z(2)-(betta/delta).*(1+pi1.*(((z(2)-z(1))./(z(1)))))];
 
Ztds = fsolve(tds,[0.1 0.1]);
 
x=Ztds(1);  % long run MRS state 1
y=Ztds(2);  % long run MRS state 2
mu=((y-x)/(x));  % long run mu
  
x_target=x   % targeted share of old consumption in ss in state 1
y_target=y   % targeted share of old consumption in ss in state 2
 
time=10; % number of iteration before convergence

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%% convergence %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% use of shooting algorithm
guess=1; % initial guess of mu0  
  
options=optimset('Display','iter','LargeScale','off','MaxIter',1000,'TolFun',1e-10,'TolX',1e-10); 
[k_sol,fval,exitflag,output,jacobian_0]=fsolve('Saddle',guess,options);
        
k=k_sol; % endogenous lagrangian multipliers mu0

% control for non-negativity constraint
if betta/delta>yaut
    y(1)=betta/delta;
else
    y(1)=yaut;
end

x(1)= ((betta/delta)/(1+k)); 

for t=1:time;
    
% MRS state 2
   if (betta/delta)*(1+pi1*(((y(t)-x(t))/(x(t)))))>yaut;
        y(t+1)=(betta/delta)*(1+pi1*(((y(t)-x(t))/(x(t)))));
    else
        y(t+1)=baut;
    end 
 % MRS state 1     
   x(t+1)=((varxi*(1+x(t))^((1/(betta*pi1)))*(((1+y(t+1))/(y(t+1))))...
       ^(((pi2)/(pi1))))/(1-varxi*(1+x(t))^((1/(betta*pi1)))...
       *(((1+y(t+1))/(y(t+1))))^(((pi2)/(pi1))))); 
 
 % Lagrangian multiplier mu   
 mu(t)=((y(t)-x(t))/(x(t)));
 
 % Consumption in state 1
 cy1(t)=e1/(1+x(t));
 co1(t)=e1-cy1(t);
 
 % consumption in state 2
 cy2(t)=e2/(1+y(t));
 co2(t)=e2-cy2(t);
 
end

x_lr=x(time);
y_lr=y(time);


plot(x,y,'r +','LineWidth',1)
hold on
plot(x_lr,y_lr,'b *','LineWidth',1)
hold on
[xx,yy]=meshgrid(0.91:.01:0.96,1.02:.01:1.18);

